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A  CLOSED  FORM  SOLUTION  OF  A  LONGITUDINAL  BEAM 
WITH  A  VISCOUS  BOUNDARY  CONDITION 


1.  INTRODUCTION 

The  dynamic  response  of  a  beam  with  a  viscous  boundary  condition  is  important 
because  the  designers  of  various  structures  often  use  viscous  dampers  to  reduce  force 
transmissibility  and  decrease  displacement.  The  need  for  reduced  force  transmissibility  is 
evident  since  lower  force  levels  permit  simpler  and  lighter  structural  designs.  For  this 
reason,  viscous  (shock)  absorbers  are  currently  critical  to  many  systems,  such  as 
buildings,  cars,  and  airplanes. 

The  closed  form  response  of  structures  with  fixed  and  free  boundary  conditions  has 
been  previously  analyzed  [1,2,3].  These  analyses  form  the  basis  for  many  classical  beam 
problems.  Their  self-adjoint  operators  are  discretized  using  mutually  orthogonal  modes  to 
produce  models  of  structural  dynamic  response.  These  models,  however,  admit  only 
standing  waves  into  the  response  and  do  not  consider  viscous  damping  at  the  boundary. 
Recently,  a  number  of  papers  have  appeared  that  model  different  effects  in  beams,  such  as 
compressive  axial  loads  [4,5],  elastic  foundations  [6],  and  the  coupling  between  flexural 
and  torsional  modes  to  axial  loads  [7].  These  papers,  which  use  a  variety  of  analytical 
techniques  to  solve  for  the  structural  response  of  a  beam,  do  not  consider  viscous 
dissipation  at  the  boundary. 

Although  the  structural  response  of  a  beam  with  a  viscous  damper  can  be 
determined  using  finite  element  analysis  [8,9],  this  method  does  have  a  number  of 
drawbacks.  It  is  computationally  intensive,  does  not  provide  explicit  eigenvalues  and 
eigenvectors,  and  does  not  yield  a  closed  form  solution.  Finite  element  discretizations  are 
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often  too  large  when  real  time  computations  are  required,  as  in  the  case  of  active  control 
systems.  Additionally,  the  effects  of  changing  model  parameters  is  not  always  as  apparent 
in  a  finite  element  model  as  it  is  in  a  closed  form  algebraic  solution. 

This  paper  develops  a  closed  form  series  solution  for  the  axial  wave  equation  with  a 
fixed  boundary  condition  at  one  end  and  a  viscously  damped  boundary  condition  at  the 
other.  The  addition  of  a  damper  to  the  boundary  allows  propagating  and  standing  waves  to 
exist  in  the  structure  simultaneously.  The  system  model  produces  a  differential  operator 
that  is  nonself-adjoint  and  corresponding  eigenfunctions  that  are  not  mutually  orthogonal. 
By  redefinition  of  the  problem  on  another  interval,  the  space  and  time  modes  will  decouple 
and  a  closed  form  series  solution  can  be  found. 

2.  SYSTEM  MODEL 

The  system  model  represents  an  axial  beam  fixed  at  x  =  0  and  a  viscous  damper  at  x 
=  L  (Figure  1).  A  force  is  applied  to  the  beam  at  location  x  =  xf.  The  addition  of  the 
damper  to  the  beam  will  admit  standing  and  propagating  wave  energy  simultaneously.  The 
linear  second  order  wave  equation  modeling  particle  displacement  in  the  beam  is 

9^u(x,t)  E  9^u(x,t)  _  5(x-Xf)F(t) 
dt^  p  dx^  pA 

where  u(x,t)  is  the  displacement  (m),  E  is  the  modulus  of  elasticity  (N/m2),  p  is  the  density 
of  the  beam  (kg/m3),  x  is  the  spatial  location  (m),  t  is  time  (s),  A  is  the  area  of  the  beam 
(m2),  F  is  the  applied  force  (N),  and  5  is  the  Dirac  delta  function  (1/m).  The  wave 
equation  assumes  uniform  area  and  negligible  internal  loss  in  the  beam. 

The  boundary  at  x  =  0  is  fixed  and  can  be  expressed  as 

u(0,t)  =  0  .  (2) 

The  boundary  condition  at  x  =  L  is  obtained  by  matching  the  strain  energy  of  the  beam  to 
the  viscous  dissipative  force  in  the  damper.  This  can  be  expressed  as 
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3u  9u 

AE— (L,t)  =  -c— (L,t) 
dx  dt 


(3) 


where  c  is  the  viscous  damping  coefficient  (Ns/m).  When  c  is  equal  to  zero  (or  infinity), 
the  boundary  at  x  =  L  reflects  all  the  wave  energy,  and  the  system  response  is  composed 
only  of  standing  waves.  When  c  is  equal  to  A^f^,  the  boundary  at  x  =  L  absorbs  all  the 

wave  energy,  and  the  system  response  is  composed  only  of  propagating  waves.  All  other 
values  of  c  exhibit  some  combination  of  standing  and  propagating  wave  energy  in  their 
response. 


3.  SEPARATION  OF  VARIABLES 

A  decoupled  series  of  ordinary  differential  equations  that  represent  the  system  are 
now  developed  from  equations  (1)  -  (3)  and  the  initial  conditions  of  the  beam.  The  first 
step  in  deriving  decoupled  differential  equations  involves  finding  the  eigenvalues  and 
eigenfunctions  of  the  model.  This  is  accomplished  by  application  of  separation  of  variables 
to  the  homogeneous  version  of  the  wave  equation  (equation  (1))  and  the  corresponding 
boundary  conditions  (equations  (2)  and  (3)).  Separation  of  variables  assumes  that  the 
independent  variable  can  be  written  as  a  product  of  two  functions:  one  in  the  time  domain 
and  one  in  the  spatial  domain.  This  form  is 

u(x,t)  =  T(t)X(x)  .  (4) 
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Equation  (4)  is  now  substituted  into  the  homogeneous  version  of  equation  (1),  which 
produces  the  following  ordinary  differential  equations; 


^-^-X^-T(t)  =  0 

dt^  P 


(5) 


and 


ff-A.^X(x)  =  0  ,  (6) 

dx^ 


where  X  is  the  complex-valued  separation  constant. 
The  general  solution  to  equation  (5)  is 


X 

T(t)  =  Ae 


(7) 


The  fixed  boundary  condition  (equation  (2)),  is  now  applied  to  the  spatial  ordinary 
differential  (equation  (5))  which  gives 

X(x)  =  e^’'-e"^’‘  .  (8) 

Applying  the  viscous  boundary  condition  (equation  (3))  to  equations  (7)  and  (8)  yields  B  = 
0  and  the  n -mode-indexed  separation  constant 


X„  = 


■~loge 

2L 


'E 

AE-cJ- 


AE-i-c,  - 

VP 


(2n  4-  l)7t  ■ 
2L  ' 


n  =0,±l,db2,.. 


(9) 


where  i  is  the  square  root  of  -1.  The  real  part  of  the  separation  constant  is  associated  with 
the  propagating  wave  energy,  and  the  imaginary  part  of  the  separation  constant  is 
associated  with  the  standing  wave  energy.  The  eigenvalues  of  the  system  are  the  separation 
constants  multiplied  by  the  wave  speed. 
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where  An  has  units  of  rad/s.  A  plot  of  the  eigenvalue  location  in  the  complex  plane  is 
shown  in  Figure  2.  The  eigenvalues  are  equally  spaced  and  parallel  to  the  real  axis.  Once 
the  indexed  separation  constant  is  determined,  the  spatial  eigenfunctions  can  be  defined  by 
inserting  equation  (9)  into  equation  (8),  which  gives 

(Pn(x)  =  e^n*-e“^"’'  .  (11) 

A  typical  eigenfunction  (n  =  2)  is  shown  in  Figure  3  for  c  =  (0.5) The 
eigenfunctions  are  not  mutually  orthogonal  on  [0,L];  therefore  their  inner  product  with 
respect  to  one  another  on  [0,L]  is  not  zero,  and  traditional  boundary  value  techniques  will 
not  decouple  the  time  and  space  modes.  A  method  is  developed  below  that  redefines  the 
problem  interval  over  [-L,L]  and  decouples  the  time  and  space  modes  of  the  system.  Once 
the  modes  have  been  decoupled,  the  problem  solution  can  be  transformed  back  to  the 
original  interval  of  [0,L]. 
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4.  SERIES  SOLUTION 

The  displacement  (or  solution)  to  the  forced  wave  equation,  written  as  a  series 
solution,  is 

oo 

u(x,t)=  ^bn(t)(pn(x)  ,  (12) 

n=-oo 

where  the  bn(t)'s  are  the  generalized  coordinates  and  the  (pn(x)’s  are  the  spatial 
eigenfunctions.  Derivation  of  an  inner  product  that  decouples  the  time  and  space  modes 
requires  the  time  derivative  (velocity)  of  the  panicle  displacement  to  be  written  in  two 
different  forms.  The  first  form  is  the  derivative  of  equation  (12)  and  yields 

-^(x,t)=  ^bn(t)(pn(x)  ,  (13) 

n=-oo 

where  the  dot  over  b  denotes  the  time  derivative  of  the  function, 
using  equations  (7)  and  (4),  is  written  as 

-^(x,t)=  2^Anbn(t)(Pn(x)  . 

n=— oo 

Equating  equations  (13)  and  (14)  produces 

oo 

^[bn(t)-Anbn(t)](Pn(x)  =  0  .  (15) 

n=— oo 

The  assumption  is  now  made  that  differentiation  will  distribute  over  the  summation. 
Decoupled  space  and  time  modes  will  validate  this  assumption.  The  forced  wave  equation 
(1)  is  rewritten  with  the  above  equations.  The  second  partial  time  derivative  is  found  from 
the  time  derivative  of  equation  (14),  and  the  second  partial  spatial  derivative  is  found  from 
the  second  spatial  derivative  of  equation  (13).  Insening  these  derivatives  into  equation  (1) 
yields 


The  second,  developed  by 


(14) 
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X  [bn(t)-Anbn(t)]An9n(x)  =  ^^^^-^^f^ 

pA 


(16) 


n=-oo 


Equation  (15)  is  now  differentiated  with  respect  to  x  and  multiplied  by  the  wave 


speed  .  The  result  is  then  added  to  equation  (16)  to  form 


X  [bn(t)-A„bn(t)]2Ane^n^=^^^^-^^i^  x  e  [0,L]  (17) 

n=—  PA 

and  is  subtracted  from  equation  (16)  to  give 

X  [bi,(t)- A„b„(t)]2Ane~’^°’‘  =  -8fa-.M)F(l)  xs(0,L]  .  (18) 

pA 

The  interval  of  equation  (18),  now  changed  from  [0,L]  to  [-L,0]  by  substitution  of -x  for  x, 
yields 


X  [bn(t)-A„b„(i)i2A„e^"’‘=  xe[-L,0J  .  (19) 

n.—  PA 

Combining  equations  (17)  and  (19)  into  a  single  equation  and  breaking  the  exponential  into 

terms  that  contain  the  index  n  and  terms  that  do  not  contain  the  index  n  results  in 

^  (2n+l)7txi 

X[bn(t)-A„bn(t)]2Ane  2L  = 

n=— oo 


—  loge 
2L 


AE-^p-cVe 

AE^Jp+c^fE 


2L 


loge 


AE-y/p-ci/E 

AE-^+cVE  J 


-5(-x-Xf)F(t) 

pA 


5(x-Xf)F(t) 

pA 


X  e  [-L,0] 


X  e  [0,L] 


(20) 
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-(2m+l)7rxi 

The  exponential  e  (where  m  is  an  integer)  is  now  multiplied  on  both  sides  of 

equation  (20),  and  the  resulting  equation  is  integrated  from  -L  to  L.  The  left-hand  side  of 

the  equation  is  orthogonal  on  the  new'  interval  and  can  be  expressed  as 

^  (2n+l)jtxi  -(2in+l)nxi 

J_Jbnm(t)-Anbn^(t)]2A„e  2L  g  2L  dx 


[bn(t)- Anbn(t)]4LAn  n  =  m 

0  n  m 


(21) 


Use  of  the  reflection  property  of  integrals  and  the  bound  of  0  <  xf  <  L  results  in  the  right 
hand  side  of  equation  (20)  becoming 


F(t) 

pA 


-5(-x-Xf)e  +  6tx-Xf)e  ^'"’‘dx 


-F(t)(pm(xf) 

pA 


(22) 


Equations  (21)  and  (22)  can  be  equated  (for  n  =  m)  to  form  ordinary  differential  equations 
for  the  generalized  coordinates  bjj  as 


bn(t)- Anbn(t)  = 


-F(t)(pn(xf) 

4LAnpA 


(23) 


An  explicit  solution  to  equation  (23)  cannot  be  found  until  a  time-dependent  forcing 
function  has  been  specified. 

The  initial  conditions  of  the  generalized  coordinates  can  be  determined  from  the 
initial  conditions  of  the  beam.  This  equation  is 
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a„(01  = 


Jo  dx  dx  4A,nLJo  3t 


(x,0)(pn(x)dx 


(24) 


where  — (x,0)  is  the  initial  strain  energy  in  the  beam  (dimensionless)  and  — (x,0)  is  the 
dx  dt 

initial  velocity  of  the  beam  (m/s).  The  formulation  of  equation  (24)  is  presented  in  the 


appendix. 


5.  FREQUENCY  RESPONSE 

A  frequency  domain  solution  to  equation  (23)  can  be  found  by  specifying  that  the 
forcing  function  F(t)  be  equal  to  a  harmonic  function  where  Fq  has  units  of 

newtons.  Solving  the  differential  equation  (23)  results  in  the  following  solution  to  the 
generalized  coordinates: 


u  _ _ Po9n(^f)  tox 

"  (iC0-An)4LAnpA 


(25) 


When  equation  (25)  is  inserted  into  equation  (12),  the  displacement  of  the  beam  becomes 


u(x,,)=  y 

"  (ico- An)  4LAnpA 


(26) 


Equation  (26)  can  be  truncated  using  2N  symmetric  terms  to  yield  an  engineering  solution 
to  the  problem. 

6.  A  NUMERICAL  EXAMPLE 

The  accuracy  of  the  model  and  the  effects  of  truncation  error  were  investigated  with 
a  numerical  example.  The  following  constants  were  used:  L  =  20  m,  E  =  207  x  10^ 
N/m^,  p  =  7.8  X  10^  kg/in^,  Fq  =  1000  N,  A  =  0.01  m^,  xf  =  3  m,  and  c  =  75000  Ns/m. 
Figure  4  is  the  frequency  domain  response  of  the  structure  u(Xi-,co)  viewed  at  Xj-  =  1 1  m  for 
a  6-term  (-3<n<2)  and  a  50-term  (-25  <  n  <  24)  model.  Numerical  simulations 
suggest  that  two  first  order  terms  are  needed  to  model  each  beam  resonance.  There  is  only 


11 


TR  10,017 


a  0.45  percent  (-46.9  dB)  difference  between  the  two  truncated  models  up  to  the  third 
resonance  of  the  beam.  The  addition  of  terms  to  the  model  does  not  change  its  value  at  the 
lower  frequencies.  This  is  due  to  the  frequency  content  of  each  generalized  coordinate:  the 
lower  indexed  coordinates  contain  only  lower  frequency  response  and  the  higher  indexed 
coordinates  contain  only  higher  frequency  response. 

Figure  5  shows  the  6-term  frequency  response  of  the  structure  compared  with  finite 
element  analysis  results  using  5-node  (square  marker),  8-node  (diamond  marker)  and  21- 
node  (triangle  marker)  finite  element  models.  The  addition  of  terms  to  the  finite  element 
model  produces  more  accurate  results  at  lower  frequencies  due  to  the  coupling  between  the 
nodes  in  a  finite  element  formulation.  Because  the  bandwidth  of  the  system  matrices  are 
greater  than  one,  the  mode  shapes  are  not  explicit  to  the  analysis,  and  the  addition  of  terms 
(nodes)  to  the  analysis  can  affect  the  accuracy  of  many  beam  resonance  modes.  The 
continuous  formulation  presented  above  eliminates  the  problem  of  banded  system  matrices 
by  use  of  an  onhogonal  inner  product  to  decouple  the  mode  shapes  of  the  beam. 
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7.  CONCLUSIONS 

The  axial  response  of  a  beam  with  a  damper  at  the  end  can  be  determined  by  using 
separation  of  variables  and  by  changing  the  interval  of  the  problem.  A  truncated  series 
solution  can  be  implemented  to  approximate  the  exact  dynamic  response.  Truncation  error 
at  lower  frequencies  is  extremely  small.  This  closed  form  solution  is  computationally 
efficient  and  the  eigenvalues  and  eigenvectors  of  the  system  are  explicit. 
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APPENDIX 

INITIAL  CONDITIONS  OF  GENERALIZED  COORDINATES 

The  initial  conditions  of  the  generalized  coordinates  can  be  determined  from  the 
initial  conditions  of  the  beam.  Rewriting  equation  (12)  at  time  t  =  0  yields 

OQ 

u(x,0)=  ^bn(0)(Pn(x)  .  (A.1) 

n=— oo 

Equation  (A.1)  is  now  differentiated  with  respect  to  the  spatial  variable  x  to  give 

^(x,0)=  ^b„(0)X„(e^n’‘+e'^"’‘)  (A.2) 

n=-oo 

and  differentiated  with  respect  to  the  time  variable  t  to  give 

|^(x,0)=  ]£b„(0)Xn(e^n’'-e'^«*)  .  (A.3) 

n=-oo 

Equation  (A.2)  is  the  initial  strain  energy  in  the  beam  and  equation  (A.3)  is  the  initial 
velocity  of  the  beam.  Equations  (A.2)  and  (A.3)  are  now  added  to  yield 

oo 

n=— oo 

and  subtracted  to  yield 

oo 

n=— oo 

The  interval  of  equation  (A.5)  is  now  changed  from  [0,L]  to  [-L,0]  by  the  substitution  of  -x 
for  X.  This  equation  is 


^(x,0)-^(x,0) 
d\  dt 


(A.5) 


^(x,0)  +  ^(x,0) 
dx  dt 


(A.  4) 


A-1 
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oo  _ 

Y  2b„(0)>Lne^"’‘ =  |^(-x.0)-|^(-x.0) 
[ox  dt 


xe[-L,0]  . 


(A.  6) 


Equations  (A.4)  and  (A.6)  are  now  combined  and  the  exponential  is  broken  into  two  terms: 
one  contains  the  index  n  and  one  does  not  contain  the  index  n.  The  resulting  equation  is 

oo  (2n+l)rtxi 

^2b„(0)>.„e  2L  = 


-1  r  AE-\/p-cVe 


rau,  AEVp+cVe 

— (-X.0)-— (-X,0)  e  L  ->  xe[-L, 

Lax  at 


(A.  7) 


-1  AE^^-c^/E 

‘au  ^  au  1 2L  ae^+cVe 
— (x,0)  +— (x,0)  e  L 

dx  at 


X  e[L,0]  . 


-(2m+l)rtxi 

The  exponential  e  jj  ^ow  multiplied  by  both  sides  and  the  left-half  of  equation 

(A.7)  is  integrated  from  -L  to  L.  This  yields 


bn(0)4XnL  =  dx 

+  f^[|^(x,0)-»-|i(x,0)le-^n’‘ dx  . 

Jo  L^x  at 


(A.8) 


Using  the  reflection  property  of  integrals  and  equation  (11),  equation  (A.8)  can  be  rewritten 


to  give 


1  _ !_ 


bn(0)  =  -^[  T^(x, 

4X,iLJo  ax 


0)-Zii^dx - - 

dx  4X 


L-f  ^ 

nLJo  at 


(x,0)(Pn(x)dx  .  (A.  9) 
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